#include "antisym_coeff.h"
int Antisym::phase(int n) { return n % 2 ? -1 : 1; }

int Antisym::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

int Antisym::Tri_state(State_Jacobi &state) {
  int unit = 1;
  int l12_2 = 2 * state.l12;
  int l3_2 = 2 * state.l3;
  if (Tri(l12_2, state.S12_2, state.J12_2) != 1) {
    return -1;
  }
  if (Tri(l3_2, unit, state.j3_2) != 1) {
    return -1;
  }
  if (Tri(state.J12_2, state.j3_2, state.J_2) != 1) {
    return -1;
  }
  if (Tri(state.T12_2, unit, state.T_2) != 1) {
    return -1;
  }
  return 1;
}

double Antisym::mat_cal(State_Jacobi &final_state, State_Jacobi &init_state) {
  double result = 0.0;
  if ((Tri_state(final_state) != 1) || (Tri_state(init_state) != 1)) {
    cout << "can`t contain Tri condition in Antisym.mat_cal";
    return 0.0;
  }

  if (phase(final_state.l12 + (final_state.S12_2) / 2 +
            (final_state.T12_2) / 2) == 1) {
    cout << "!!!WRONG!!!" << endl;
    exit(0);
  }
  if (phase(init_state.l12 + (init_state.S12_2) / 2 + (init_state.T12_2) / 2) ==
      1) {
    cout << "!!!WRONG!!!" << endl;
    exit(0);
  }

  if ((final_state.J_2 != init_state.J_2) ||
      (final_state.T_2 != init_state.T_2))
    return result;

  int E12_final = 2 * final_state.n12 + final_state.l12;
  int E3_final = 2 * final_state.n3 + final_state.l3;

  int E12_init = 2 * init_state.n12 + init_state.l12;
  int E3_init = 2 * init_state.n3 + init_state.l3;

  if ((E12_final + E3_final) != (E12_init + E3_init))
    return result;
  int parity_final = final_state.l12 + final_state.l3;
  int parity_init = init_state.l12 + init_state.l3;
  if (phase(parity_final) != phase(parity_init))
    return result;

  int phase_final =
      (final_state.S12_2 + final_state.T12_2) / 2 + final_state.l12;
  if (phase(phase_final) != -1)
    return result;
  int phase_init = (init_state.S12_2 + init_state.T12_2) / 2 + init_state.l12;
  if (phase(phase_init) != -1)
    return result;

  int L_min = std::max(std::abs(final_state.l12 - final_state.l3),
                       std::abs(init_state.l12 - init_state.l3));
  int L_max = std::min((final_state.l12 + final_state.l3),
                       (init_state.l12 + init_state.l3));
  if (L_min > L_max)
    return result;

  int S2_min =
      std::max(std::abs(final_state.S12_2 - 1), std::abs(init_state.S12_2 - 1));
  int S2_max = std::min((final_state.S12_2 + 1), (init_state.S12_2 + 1));
  if (S2_min > S2_max)
    return result;

  // part1 ~~~
  double part1 =
      sixJ.quick(1, 1, final_state.T12_2, 1, final_state.T_2, init_state.T12_2);
  if (part1 == 0)
    return result;

  // part2 ~~~
  double part2_t = (final_state.T12_2 + 1) * (final_state.J12_2 + 1) *
                   (final_state.j3_2 + 1) * (final_state.S12_2 + 1) *
                   (init_state.T12_2 + 1) * (init_state.J12_2 + 1) *
                   (init_state.j3_2 + 1) * (init_state.S12_2 + 1);
  double part2 = std::sqrt(part2_t * 1.0);

  // part3 ~~~
  double d_hobar = 1.0 / 3.0;
  double part3 = 0.0;
  for (int S2 = S2_min; S2 <= S2_max; S2 = S2 + 2) {
    for (int L = L_min; L <= L_max; L++) {
      double part3_1 = (S2 + 1) * sixJ.quick(1, 1, final_state.S12_2, 1, S2,
                                             init_state.S12_2);
      if (part3_1 == 0)
        continue;
      // double part3_2=0.0;

      double part3_2t =
          (2 * L + 1) * ninJ.cal(2 * final_state.l12, final_state.S12_2,
                                 final_state.J12_2, 2 * final_state.l3, 1,
                                 final_state.j3_2, 2 * L, S2, final_state.J_2);
      if (part3_2t == 0)
        continue;
      part3_2t =
          part3_2t * ninJ.cal(2 * init_state.l12, init_state.S12_2,
                              init_state.J12_2, 2 * init_state.l3, 1,
                              init_state.j3_2, 2 * L, S2, init_state.J_2);
      if (part3_2t == 0)
        continue;
      part3_2t =
          part3_2t * hobra_x.cal(E12_final, final_state.l12, E3_final,
                                 final_state.l3, E12_init, init_state.l12,
                                 E3_init, init_state.l3, L, d_hobar);
      // double d=1/3.0;
      // double hob=phase(final_state.l12+final_state.l3)		\
	//*hobra_x.cal(E3_final,final_state.l3,E12_final,final_state.l12, \
	//	     E3_init, init_state.l3,E12_init, init_state.l12, L,d);
      // part3_2t = part3_2t * hob;
      // if(part3_2t == 0)continue;
      // cout<<"###  "<<"L = "<<L<<"\t"<<part3_2t<<endl;
      part3 += part3_1 * part3_2t;
    }
  }

  int phase_t = phase(final_state.l12 + init_state.l12);
  // phase_t=1;
  if (-1 == 0) {
    cout << "@@@!!!@@@" << endl;
    cout << phase_t << "\t" << part1 << "\t" << part2 << "\t" << part3 << endl;
  }
  result = phase_t * part1 * part2 * part3;
  return result;
}

void Antisym::JTmat_build(Jacobi &Jacobi_t, int channel_num) {
  int size = Jacobi_t.state[channel_num].size();

  this->anti_mat.resize(size, size);

  // cout<<"!!! size = "<<size<<endl;
  //#pragma omp parallel for
  int i_x = -1;
  for (int i = 0; i < size; i++) {
    // cout<<i<<endl;
    for (int j = 0; j < size; j++) {
      // cout<<i<<"\t"<<j<<endl;
      anti_mat(i, j) = -2 * mat_cal(Jacobi_t.state[channel_num][i],
                                    Jacobi_t.state[channel_num][j]);
      //
    }
    // this->anti_mat(i,i)=1.0/3.0+this->anti_mat(i,i);
  }

  for (int i = 0; i < size; i++) {
    anti_mat(i, i) = 1 + anti_mat(i, i);
    for (int j = 0; j < size; j++) {
      anti_mat(i, j) = 1.0 / 3.0 * anti_mat(i, j);
    }
  }
}
void Antisym::JTmat_build(Jacobi &Jacobi_t, int channel_num, int E) {
  int size_C = Jacobi_t.state[channel_num].size();
  int size_E = 0;
  for (int i = 0; i < size_C; i++) {
    if (E == Jacobi_t.state[channel_num][i].E)
      size_E++;
  }

  this->anti_mat.resize(size_E, size_E);

  // cout<<"!!! size = "<<size<<endl;
  //#pragma omp parallel for
  int i_x = -1;
  for (int i = 0; i < size_C; i++) {
    // cout<<i<<endl;
    if (Jacobi_t.state[channel_num][i].E != E)
      continue;
    i_x++;
    int j_x = -1;
    for (int j = 0; j < size_C; j++) {
      // cout<<i<<"\t"<<j<<endl;
      if (Jacobi_t.state[channel_num][j].E != E)
        continue;
      j_x++;
      anti_mat(i_x, j_x) = -2 * mat_cal(Jacobi_t.state[channel_num][i],
                                        Jacobi_t.state[channel_num][j]);
      //
    }
    // this->anti_mat(i,i)=1.0/3.0+this->anti_mat(i,i);
  }

  for (int i = 0; i < size_E; i++) {
    anti_mat(i, i) = 1 + anti_mat(i, i);
    for (int j = 0; j < size_E; j++) {
      anti_mat(i, j) = 1.0 / 3.0 * anti_mat(i, j);
    }
  }
}

void Antisym::JTmat_build(Jacobi &Jacobi_t, int channel_num, int E_min,
                          int E_max) {
  int size_C = Jacobi_t.state[channel_num].size();
  int size_E = 0;
  for (int i = 0; i < size_C; i++) {
    int E_x = Jacobi_t.state[channel_num][i].E;
    if (E_x >= E_min && E_x <= E_max)
      size_E++;
  }
  if (cout_flag == 1) {
    cout << " channel " << channel_num << "  E" << E_min << endl;
    cout << "!!! size_E = " << size_E << endl;
  }
  if (size_E > 0) {
    build_flag = 1;
    this->anti_mat.resize(size_E, size_E);
  } else {
    build_flag = -1;
    this->anti_mat.resize(0, 0);
    return;
  }

  //#pragma omp parallel for
  int i_x = -1;
  for (int i = 0; i < size_C; i++) {
    // cout<<i<<endl;
    int E_x = Jacobi_t.state[channel_num][i].E;
    if (E_x < E_min || E_x > E_max)
      continue;
    i_x++;
    int j_x = -1;
    for (int j = 0; j < size_C; j++) {
      // for(int j=0;j<(size_C/2+1);j++){
      // cout<<i<<"\t"<<j<<endl;
      int E_y = Jacobi_t.state[channel_num][j].E;
      if (E_y < E_min || E_y > E_max)
        continue;
      j_x++;
      if (j_x <= i_x + 1) {
        anti_mat(i_x, j_x) = -2 * mat_cal(Jacobi_t.state[channel_num][i],
                                          Jacobi_t.state[channel_num][j]);
      } else {
        anti_mat(j_x, i_x) = anti_mat(i_x, j_x);
      }
      //
    }
    // this->anti_mat(i,i)=1.0/3.0+this->anti_mat(i,i);
  }

  for (int i = 0; i < size_E; i++) {
    anti_mat(i, i) = 1 + anti_mat(i, i);
    for (int j = 0; j < size_E; j++) {
      anti_mat(i, j) = 1.0 / 3.0 * anti_mat(i, j);
    }
  }
}

void Antisym::mat_diag() {
  // SelfAdjointEigenSolver<MatrixXd> anti_solve(anti_mat.topLeftCorner(20,20));
  if (build_flag == -1) {
    if (cout_flag == 1)
      cout << "bulid error" << endl;
    diag_flag = -1;
    anti_basis_num_d = 0;
    return;
  } else {
    int size = anti_mat.cols();
    // for (int i = 0; i < size; i++) {
    //   for (int j = 0; j < size; j++) {
    //     cout << i << " \t " << j << "\t" << anti_mat(i, j) << endl;
    //   }
    // }
    SelfAdjointEigenSolver<MatrixXd> anti_solve(anti_mat);

    eig_val_t.resize(size);
    for (int i = 0; i < size; i++) {
      // cout<<i<<"\t"<<anti_solve.eigenvalues()(i)<<endl;
      eig_val_t[i] = anti_solve.eigenvalues()(i);
    }

    anti_basis_num_d = (anti_solve.eigenvalues().sum());
    anti_basis_num = std::floor(anti_solve.eigenvalues().sum() + 0.001);
    if (std::abs(anti_basis_num_d - anti_basis_num) > 0.01 ||
        anti_basis_num_d < 0.01) {
      if (cout_flag == 1)
        cout << " Partial antisymmetry Jacobi Basis is incompleteness!!! "
                "anti_basis_num_d = "
             << anti_basis_num_d << endl;
      // exit(0);
      diag_flag = -1;
      anti_vec_C = anti_solve.eigenvectors().rightCols(anti_basis_num);
      anti_vec_D = anti_vec_C * anti_vec_C.transpose();
    } else {
      // cout<<"!!!!"<<anti_basis_num<<endl;
      anti_vec_C = anti_solve.eigenvectors().rightCols(anti_basis_num);
      anti_vec_D = anti_vec_C * anti_vec_C.transpose();
      diag_flag = 1;
    }
    return;
  }
}

void Antisym::cal() {
  // cout<<"&&^^&&"<<endl;
  // cout<<Jacobi_x.channel.size()<<endl;
  this->cout_flag = -1; // 1: print details  others don't print

  this->eig_val.resize(Jacobi_x.channel.size());
  this->vec_D.resize(Jacobi_x.channel.size());
  this->eig_num.resize(Jacobi_x.channel.size());

  int complete_flag = 1; // flag for completeness of Jacobi basis

  double plan_t = 0.0;

  int pos = 0;
  if (cout_flag == 1) {
    cout << "Begin to calculate antisymmetry Jacobi basis" << endl;
  }

  Pgs_Bar pgs;
  pgs.init(Jacobi_x.channel.size(), 30);
  //#pragma omp parallel for schedule(dynamic, 1)
  for (int i = 0; i < Jacobi_x.channel.size(); i++) {
    if (cout_flag == 1) {
      pgs.PrintBar(i);
    }

    eig_val[i].resize(Jacobi_x.channel[i].E_max + 1);
    vec_D[i].resize(Jacobi_x.channel[i].E_max + 1);
    eig_num[i].resize(Jacobi_x.channel[i].E_max + 1);

    complete_flag = 1;
    for (int E = 0; E <= Jacobi_x.channel[i].E_max; E++) {

      int E_temp = std::max(1, E);
      if ((E > Jacobi_x.channel[i].E_min + 1) && (eig_num[i][E_temp - 1] < 0)) {
        complete_flag = -1;
      } else {
        if (cout_flag == 1)
          cout << " channel : " << i << " E : " << E
               << "\t c_n : " << Jacobi_x.channel[i].num << endl;
        this->JTmat_build(Jacobi_x, i, E, E);
        this->mat_diag();
      }

      if (build_flag == -1 || diag_flag == -1 || complete_flag == -1) {
        this->eig_val[i][E].resize(0);
        this->vec_D[i][E].resize(0);
        this->eig_num[i][E] = -1.0;
        // cout<<"bad !! bad"<<endl;
      } else {
        // cout<<"E = "<<E<<endl;
        this->mark_Jaco(i, E);
        int size_E = anti_vec_D.cols();
        this->eig_val[i][E].resize(size_E);
        this->vec_D[i][E].resize(size_E);
        eig_num[i][E] = anti_basis_num_d;
        // cout<<"size_E : "<<size_E<<"\t anti_vec_D size
        // "<<anti_vec_D.size()<<endl;
        for (int x = 0; x < size_E; x++) {
          eig_val[i][E][x] = eig_val_t[x];
          vec_D[i][E][x].resize(size_E);
          for (int y = 0; y < size_E; y++) {
            // vec_C[i][E][x][y]=anti_vec_C(x,y);
            vec_D[i][E][x][y] = anti_vec_D(x, y);
          }
        }
      }
    }
  }
  if (cout_flag == 1) {
    pgs.EndBar();
    cout << "Finish to calculate antisymmetry Jacobi basis" << endl;
  }
}

void Antisym::mark_Jaco(int channel, int E) {
  // mark the flag of Jacobi basis to 1 for completeness; -1 for uncompleteness
  // (default)
  for (int i = 0; i < Jacobi_x.state[channel].size(); i++) {
    if (Jacobi_x.state[channel][i].E == E) {
      Jacobi_x.state[channel][i].flag = 1;
    }
  }
}

void Antisym::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = vec_D.size();
  for (int i = 0; i < size_ch; i++) {
    file << "channel : " << i << endl;
    int size_E = vec_D[i].size();
    for (int E = 0; E < size_E; E++) {
      file << "=== E : " << E << "  ===" << endl;
      int size = vec_D[i][E].size();
      for (int j = 0; j < size; j++) {
        for (int k = 0; k < size; k++) {
          file << j << "\t" << k << "\t" << vec_D[i][E][j][k] << endl;
        }
      }
    }
  }
  file.close();
}
